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We discuss the Meissner response to a known field source of superconductors having inhomo- 
geneities in their penetration depth. We simplify the general problem by assuming that the pertur- 
bations of the fields by the penetration depth inhomogeneities are small. We present expressions for 
inhomogeneities in several geometries, but concentrate for comparison with experiment on planar 
defects, perpendicular to the sample surfaces, with superfiuid densities different from the rest of 
the samples. These calculations are relevant for magnetic microscopies, such as Scanning Supercon- 
ducting Quantum Interference Device (SQUID) and Magnetic Force Microscope, which image the 
local diamagnetic susceptibility of a sample. 
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I. INTRODUCTION 

Recently Kalisky et al. reported the observation of 
"stripes" of enhanced local diamagnetic susceptibility 
in scanning SQUID microscope (SSM) images of un- 
derdoped Ba(Fei_; C Co :r )2As2pJ They interpreted these 
stripes as being due to enhanced superfiuid density (re- 
duced London penetration depth) along twin boundaries. 
Kirtley et al. modeled these experiments by solving Lon- 
don's and Maxwell's equations using finite element meth- 
ods in an appropriate geometryP These calculations were 
in agreement with the experimental lineshapes, and pro- 
vided estimates for the increase in the superfiuid density 
on the twin planes. However, it was difficult to model ge- 
ometries with regions of enhanced superfiuid density as 
narrow as seemed physically likely, so that extrapolations 
to narrow widths from wider ones were necessary. 

The geometry relevant to this problem, that of a slab 
imbedded perpendicular to the sample surface in a bulk 
half-space superconductor, is difficult to treat analyti- 
cally. Here we make the problem tractable by (1) assum- 
ing that the width of the region with reduced penetration 
depth is small relative to other lengths in the problem, 
and (2) treating the problem to first order in a pertur- 
bation expansion. The first assumption is most likely 
valid for the case of SSM, since in this case the experi- 
mentally observed stripes in susceptibility are resolution 
limited.1' Treating the problem to first order in a pertur- 
bation expansion seems reasonable, since at least at low 
temperatures the stripes in susceptibility observed using 
SSM are much smaller than the susceptibility itself. 



A. SSM technique 

Although the method developed here for the evaluation 
of the Meissner response of superconductors with inho- 
mogeneities is general, we will use as a concrete example 
scanning SQUID susceptometrypl which employs a sensor 
with two concentric, co-planar loops: one loop carries a 



small current / that is a source of a weak magnetic field, 
and the other loop couples the response magnetic flux 
into the sensor SQUID. This is an elaboration of the com- 
mon SQUID magnetometry, in which a SQUID senses the 
intrinsic magnetic fields without a source coil. Our re- 
sults are relevant for magnetic force microscopy (MFM) 
as well. However, applying our approach to MFM re- 
quires modeling of magnetic tips with complex geometry 
and is outside of the scope of this paper. 

Twin and grain boundaries in superconductors may 
have enhanced as well as suppressed superfiuid density. 
We will use below a generic term "defect". The theory 
developed here applies for both enhancement and sup- 
pression, provided that the deviation of the superfiuid 
density (or of the London penetration depth) at the de- 
fect from the bulk value is small. 



B. Method 

Let us consider a magnetic field source with known field 
distribution h s in the absence of a superconductor. The 
source is placed above the superconducting half-space 
z < 0. The total field in the empty half-space z > 
can be written as 

h = h s + h r , (1) 

where h r is the response field, which satisfies div h r = 
curl h r = in vacuum outside the superconductor. One 
can look for this field as V</J r , with the scalar potential 
tp r obeying the Laplace equation and the boundary con- 
dition that it approaches zero far from the surface. The 
general form of such a potential is 

^ (r ' z) = / (0^ P(fc)ea, ' r ~**- (2) 

Here, r = (x, y), k — (k x , k y ), and z is directed normal to 
the superconducting flat surface at z — 0; Lp r (k)e~ kz is 
the two-dimensional (2D) Fourier transform with respect 
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to variables x, y at any fixed z > 0. The potential (2j is 
denned only in the upper half-space; hence, the problem 
of uniqueness that is in general associated with the de- 
scription of the static magnetic field by a scalar potential 
does not arise. 

Thus, to know the outside field distribution it suffices 
to find the 2D Fourier transform ip r (k). In principle, this 
can be done by solving the internal London problem and 
by utilizing the boundary conditions of field continuity at 
the interface z = 0. The formal difficulty to overcome is 
to calculate the 2D Fourier transform of internal fields for 
non-uniform superconductors. Below, we show how this 
can be done in a few cases relevant for SSM and MFM. 



C. Uniform and isotropic half-space 

To demonstrate the method, we start with the sim- 
ple situation of a uniform and isotropic half space, for 
which the London equation is Hq — AqV 2 /io = 0, with 
Ao being the London penetration depth. The 2D Fourier 
transform then reads: 



p 2 h Q (k,z) -h%(k,z) = 0, p 2 = \^ 2 + k 2 



(3) 



where the prime denotes d/dz. The solution that van- 
ishes at z — > — oo is 



ho(k,z) = H(k)e pz 



(4) 



with H independent of z. h should satisfy div h Q = 0, 
which yields in Fourier space: 

i{k x H x + k y H y ) + P H Z = . (5) 

The requirement of field continuity at z = gives: 

H x = ik x ((p s + ip%) , (6) 

Hy = iky(<p S + ifl) , (7) 

H z = k(tp* - p£) . (8) 

We took into account here that for the source of the mag- 
netic field placed at z = z , the potential under it, in 
particular at z = 0, is given by 



<P 8 (r,z) = 



d 2 k 



tp s (k)e lk - r+k{z - Zo) , 



so that h s z (k) = +k<£_ s (k) 

Multiplying Eq 
adding them up, an 



fby ik x , ((7k by ik y , and |[8) 
using Eq.l5l) yields: 



^(fc) = ^V(fc) 



p + k 



TT 2ik X: yP 



rr 2fc 2 s 

p + k 



(9) 

by p, 

(10) 
(11) 



Thus, the response fields outside and inside are expressed 
in terms of the unperturbed source field tp s . This result 
has been obtained in Ref.|4]as a particular case of cumber- 
some anisotropic formulas; here it follows directly from 
the isotropic London equations. 



II. PLANAR DEFECT 

For the general case of an inhomogeneous penetration 
depth A(r), the magnetic field within the superconductor 
obeys the London equation in the form: 

47T 47T 47T 

h H curl(A 2 j) = h H A 2 curl j + VA 2 x — j = 0. 

c * c c 

(12) 

This equation is the minimum condition for the London 
energy functional 



£ l = J ^[h 2 + X 2 (cm\hf 



(13) 



which holds for inhomogeneous A. 

For a planar defect at x = 0, we model the penetration 
depth by 



A 2 (x) = X 2 -p 3 S(x) 



(14) 



where a positive /3 with the dimension of length is related 
to a superfluid density enhancement, whereas j3 < cor- 
responds to a superfluid density suppression. Physically, 
the superfluid density at the planar defect may extend 
to distances on the order of the coherence length £ into 
the bulk. However, within the London approach for ma- 
terials with £ <C A the representation ( 14 1 is justified. 



The advantage of Eq. ( 14 ) is that it allows one to do the 



2D Fourier transform of the London equation for which 
analytic expressions for all transformed quantities on the 
whole x, y plane are needed. 

An alternative way to address the problem could be to 
consider the defect as a layer of a finite thickness with the 
penetration depth different from Ao of the surrounding 
material, to look for solutions of the London equations 
in each part separately and to match them with certain 
boundary conditions. These real space solutions should 
then be matched with the real space field distribution 
in the outer space to calculate the response field. This 
approach, however, is more cumbersome and certainly 
less tractable and transparent as compared to the method 
utilizing the 2D Fourier transform employed here. 

With \(x) of Eq. (14), the London equation ( 12 ) takes 
the form 

h - \lV 2 h = f3 3 6'(x)x x curl h - (3 3 6(x)V 2 h. (15) 

The idea of the following manipulation is based on the 
physical assumption that the influence of the defect on 
the field distribution is weak, j3 <C Ao, and one can use 
a perturbation argument for its evaluation. Fits of the 
present theory to the experiments of Kalisky et al. (Fig. 
[3]) require values of /3 ~ A . However, comparison of fi- 
nite element modeling of the same problem (Fig. [2]) are 
in reasonable agreement with the present theory, even 
for /3 ~ Ao. This justifies keeping only the first order in 
perturbation theory, resulting in a considerable simplifi- 
cation of the problem. 
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Having this in mind, we look for the held inside as 
h = Hq + h b , where the unperturbed held satishes Hq — 
AgV 2 fa = in the absence of the defect plane, whereas 
h b is a perturbation due to the boundary. Wc then obtain 
in the hrst order: 

h b - X 2 V 2 h b = f3 3 6'{x)x x curlh - (3 3 5{x)V 2 h a , (16) 

where ho has been calculated in the preceding section. 

One now calculates the 2D Fourier transform (FT) of 
the left-hand side (LHS): 

FT(h b - \ 2 V 2 h b ) = (1 + X 2 k 2 )h b {k, z) - \ 2 h' b \k, z). (17) 



use easily verifiable identities, see Appendix A 

dq x 



Calculating the 2D FT of the RHS of Eq. fl6|), one can 

, see Appi 

f\Qx> ky) 



FT[5(x)f(r)] = 
FT[6'(x)f(r)]=i 



, 2^ 
°° dq x 



(18) 



3G 



2tt 



(kx Qx 

)f(q x> k y ). (19) 



We obtain after straightforward algebra: 

P 2 h b {k,z)-K(Kz)=^ r 



(20) 



where out of the three components of the vector A we 
will need only one: 

A z = (ih' 0x + q x h 0z )(kx - q x ) - Kz + Q 2 hoz , (21) 

where Q = (q x , k y ). The held h (Q,z) satisfies Eq. Q 
in which one should replace k > Q and p — » K: 



ti>(Q,z) = K z h Q (Q,z), K=J\~ 



Hence, 



h {Q,z)=H(Q)e 



Kz 



with H given in Eq. (Ill: 



H *' v ~ K + Q 



2iQ x , y K sr ^ a 2Q J 



2 + Q 2 . (22) 
(23) 

-y(Q). (24) 



Substituting these H(Q) in Eq. pl| ) we obtain: 
A z = -2{K - Q)(Q • k)e Kz V s {Q) = A 0z e Kz 



(25) 



We now write Eq. <|20j) for the field perturbation in a 
compact form: 



h' b \k,z) -p 2 h b {k,z) = D(k,z), 
D(k,z) = -^ ^ ^e^A . 



(26) 



This is a second order linear differential equation for 
h b (k, z) with respect to the variable z. The solution van- 
ishing at z — > — oo is 



h b (k,z) = Ce p: 



A 2 , 



Kz 



J-oo 



dq x e 
^2tt K 2 - p 2 



A (27) 



(see Appendix B). The arbitrary vector C = (C x , C y ,C z ) 
is to be determined from the boundary conditions. 

In fact, the constants Cj are not independent because 
div/if, = 0. In particular, at z — this gives 



ik ■ C 



/3 3 



dq x (ik ■ A + KA Qz ) 



2tt 



K 2 



P" 



(28) 



Now, we can formulate the boundary conditions of field 
continuity at z = 0: 



ik x (ip s 

iky((p S 



¥> r ) = h Qx 

■<P r ) 
ip s ) = h a 



Oy 



c z - 



P 3 


/• X 


dq x 


Aq x 




/ — oo 


2tt 


I 2 ~ k 2 


/3 3 


p oo 


dq x 


A 0y 




' — oo 


2tt 


I 2 ~ k x 


£j 


f> X 


dq x 


A 0z 


^0 - 


— oo 


2ir 


<& - k l 



(29) 
(30) 
(31) 



Multiply the hrst equation by ik x , the second by ik y , 
and the third by p and add them up. The terms with 



h add to zero because divh = 0. Utilizing Eq. (28) we 
obtain for the defect contribution to the outside magnetic 
potential: 



_ 2f3 3 
~ \ 2 k(k+p) 



<p r (k) = V r {k) 



P 



p + k 
dq x (K - Q)k ■ Q 



2tt 



K 



P 



(32) 



III. APPLICATION TO SQUID 
SUSCEPTOMETRY 



The potential of a circular current source of the SQUID 
susceptometer is given by its 2D Fourier transform!^ 



V s (k) = 



47r 2 ia 
ck 



— kzo 



J\(ka)e 



-ik x XQ 



(33) 



where / is the current through the field coil of radius a, 
(xq, 0, Zo) are the coordinates of the coil center, and zq is 
the height of the coil above the sample surface. 



A. Uniform sample 

The potential of the response field is given in Eq. ( [To| , 
so that the 2D FT of the response field ho z (k) — — kip^k) 
for a superconducting half-space free of defects is given 
by 



h 0z (k,z ) 



At: 2 la p — k 



kz 



P 



k 



Jiika) : 



(34) 



here we have set xq = since all positions xq are equiv- 
alent in this case. This gives the distribution of the z 
component of the field in the SQUID plane: 



h 0z (r,z ) 



la fd k(p- k) eikr _ 2kzo 
p + k 
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FIG. 1. (Color online) Calculated SQUID susceptometer 
response fields: a) —cho z (x,y)/4:iT 2 Ia, where ho z is the z- 
component of the field at the SQUID plane, I is the current 
through the field coil of a radius a, in the absence of a planar 
defect, b) — cX^htz/^ivip 3 , where hbz is the z-component of 
the field due to the planar defect, Ao is the London penetra- 
tion depth of the bulk superconductor, and /3 determines the 
size of the change in penetration depth at the planar defect, 
Eq. (141, for xo/a = 0, c) xo/a — 1/2, and d) xo/a — 1. Here 
z /a = 0.17 and A/a = 0.05. 



This distribution is shown in Fig.[T£i for the parameters 
indicated in the caption. Integrating this over the SQUID 
loop area of radius r, we obtain the flux of the response 
field: 



■nlar 



dk 



P 



e- 2kZ0 Ji(jfco)Ji(A:r) ; (36) 



where p 2 = A 2 + k 2 . 



B. Planar defect 

The Fourier transform of the z-component of the re- 
sponse field due to the planar defect is given by Sh z (k) = 
— kip(k) with -0 given in Eq. (32 1 and f s {Q) obtained 
from Eq. (33 1 with k replaced by Q: 



dq x ^ - Q)(fc ' Q) Jl(Qa) c - lQ+k)za - iq ' X0 



(37) 



(k+p){K+p)Q 

Q = {q x ,k y ), p=Jk 2 - 



1/X 2 . 



Here the integration over q x is done numerically for each 
k and the results are Fourier transformed to obtain the 
magnetic fields as a function of position in real space. 
Selected results for the fields are shown in Figure [T] Fig. 
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FIG. 2. (Color online) Calculated normalized SQUID sus- 
ceptibility 5%/xoi where 8\ is the change in the SQUID sus- 
ceptibility due to the planar defect, and xo is the SQUID 
susceptibility in the absence of a planar defect, for yo = 0, 
with r — 0.25a, zq = 0.17a: The curve labelled "Current 
model" (open symbols) evaluates Eq.'s (34 1 and (371 with 
P 3 /Xla = 0.15. The curve labelled "Finite element" (closed 
symbols) is numerical modeling of a stripe with finite width 
w/a — 0.2 with Xo/a = 0.2 and At/a = 0.1. 



[T^, shows the response field ho(r,z = 0), for a bulk su- 
perconductor in the absence of a planar defect. Negative 
response fields (colored red) correspond to diamagnetic 
shielding. Fig.[T]D-d display the change in the response 
field, hbz(r,zo), due to a planar defect at various spac- 
ings Xo between the center of the field coil and the defect 
position. 

Next, hbz{T,Zo) is integrated numerically over the 
SQUID loop of a radius r centered at (xo,0) to obtain 
the change in magnetic flux The integration can also 
be done analytically, see Appendix C: 



2tt 



dk 

T 



h bz (k,z )Ji(kr)e lk - x °. 



(38) 



SQUID susceptibilities are defined as \ — ®b/I$Q, 
where <&o = h/2e is the superconducting flux quantum. 
The curve labelled "Current model" in Fig. [2] shows the 
change in susceptibility 8\ due to a planar defect at x — 
divided by the susceptibility \o i n the absence of a defect 
as a function of the position xq of the SQUID sensor, with 
fixed y = 0, z /a = 0.17, r/a = 0.25, f3/a = 0.18 and 
Ao = 0.025. The parameters Ao and j3 were chosen for 
convenience of comparison with finite clement modeling 
to be discussed in Section lill CI 

Figure [3] displays the predicted S\/xo versus xq along 
with the data of Kalisky et al. taken on a twinned crystal 
of Ba(Fei_ K Co a; )2As2!^' In this case, the fitting parame- 
ters were the positions of the twin boundaries, an overall 
scaling factor (corresponding to adjusting /3 3 /AQa), the 
field coil radius a, and a vertical shift of the data. The 
fixed parameters were zo/a = 0.17 and Xo/a = 0.05. 
The agreement between experiment and theory is rea- 
sonable. The double maxima structure predicted by the 
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ting a = —i/^oX 2 and u> = 1 and recognizing that the sec- 
ond term on the LHS of Eq. ( 39 ) is quite small. We used 
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FIG. 3. (Color online) Fit of current model to experimental 
data of Ref.Q] The fitting parameters are the positions in x 
of the 8 peaks in susceptibility, an overall vertical shift, the 
radius of the field coil a, and a vertical scaling factor f3 3 /Xqcl. 
The fixed parameters are zo/a = 0.17 and r/a = 0.25, where 
r is the radius of the SQUID pickup loop. The best fit values 
area = 7.1(4-2-1.4) fj,m and P 3 /X%a = 0.048(+0.010-0.011), 
using a doubling of the best fit chi-square value as the criterion 
for determining the uncertainty in the fit parameters. 



theory is not observed in experiment, although this may 
be due to an insufficient signal-to-noise ratio. Also, we 
model the current ring and SQUID loop by linear cir- 
cles whereas both of them have a width on the order 
of microns, making it difficult to resolve the structure 
of x( x o) ° n this length scale. The data are best fit by 
0.037 < /3 3 Ao« < °- 058 and 5.7^m < a < 9.1 ^m, us- 
ing a doubling of the best-fit chi-square as a criterion 
for judging the uncertainty in the parameters. The best 
fit value for a is consistent with the inside radius of 
6/itm and outside radius of 11.5 fim of the field coil used 
in this experiment. If we assume a penetration depth 
of A = 0.325 fim for Ba(Fei_ a; Co a ;)2As2, then 0.28 
< /3 < 0.38 ^m: ~ A. 



C. Comparison with a finite element calculation 

As a consistency check, we compare our results with 
those of a finite element calculation using the commercial 
software package Comsol with the ac/dc module. This 
module solves the equation of electromagnetism in con- 
ducting media (in S.I. units) 



(iuicr — w 2 eoe r )A + curlB//io/^r = Je 



(39) 



where a is the conductivity, eo and jjbQ are the permit- 
tivity and permeability of vacuum, e r and [i r are the 
relative permittivity and permeability, B = curl A, and 
J e is an external current source. Eq. (39) can be trans- 
formed into London's equation V 2 A — A/ A 2 = by set- 



this procedure to solve the problem of a stripe of width 
w and penetration depth At, centered at x — xq imbed- 
ded in a bulk superconductor with penetration depth Ao 
occupying the half-space z < (with fi r = 1). The field 
coil is modeled as a torus centered at [0, 0, Zq] with ma- 
jor radius a and minor radius b/a = 0.05. The boundary 
conditions were continuity of A at the internal bound- 
aries and n x A = at the external boundaries. The 
results of this calculation were qualitatively similar to 
those obtained in Ref.0 although in that work the Lon- 
don equation V 2 B — A 2 B = was solved, resulting in 
solutions that did not necessarily satisfy the condition 
divB = 0. The current finite element calculations solve 
the London equation for the vector potential A, assuring 
that divB = divcurlA = 0. 

A plot of the resultant normalized change in suscepti- 
bility 5x1 X due to a slab of superconductor with width 
w/a = 0.2, Xb/a = 0.1, imbedded in a bulk supercon- 
ductor with Ao/a = 0.2 is shown as the curve labelled 
"Finite element" (solid symbols) in Fig. [2] Here the pa- 
rameters zo/a = 0.17, and r/a = 0.25 were used. For 
comparison we scaled the current model predictions using 
/3 3 = (A 2 , - X 2 b )w = 0.182a, so that ^/X^a = 0.15, (open 
symbols) in Fig. [2] The finite element calculation pro- 
duces a broader lineshape than the current model. This 
is to be expected because of the finite width assumed for 
the region of reduced penetration depth. Aside from this 
difference the two sets of results are comparable, even 
though /? 3 /Aq = 0.75 is not much smaller than 1. These 
results justify, at least in retrospect, our use of the first 
order perturbation theory in the current work. Note that 
the minimum at Xq = shown in Fig. [2] in both the cur- 
rent model and in Comsol calculation has not been seen 
experimentally, presumably due to insufficient resolution. 
We thus postpone discussion of this minimum until im- 
provements in experimental techniques make it relevant. 

We note in this connection that in our model we keep 
terms of the order /3 3 /Aq and neglect the terms C(/3 6 /Ag) 
(take, e.g., Eq. ( p2| , use Ao as a unit of length to make 
the integral dimcnsionless and on the order of 1 to see 
that V ~ /3 3 /Ag). 



IV. POINT DEFECT AND GREEN'S 
FUNCTION 

Consider a defect as a "vertical" line crossing the inter- 
face at a point r and extending from z = to z = - oo. 
Physically, such a defect affects the outside response only 
from the depth on the order of the penetration depth A. 
The rest of the defect line is irrelevant, so that one can 
consider the defect line as uniform along z. This is, of 
course, a restriction, but it allows us to treat the pene- 
tration depth as two-dimensional and to model it as 



A 2 (r) = A 2 -r7 4 <5(r-r ). 



(40) 



G 



where r\ is a constant with the dimension of length. The 
solution then provides a Green's function for a general 
problem of arbitrary distribution of such defects close to 
the sample surface relevant for the SSM technique. 
We have instead of Eq. ( 15 |: 



h b - X 2 W 2 h b = -r, A 5{r - r )W 2 h 

+ i] i S(y — y )5 / (x — x )x x curlh, 

+ r] 4 S(x - x )S'(y - y )y x curl/i . (41) 

Here, primes denote derivatives of the delta-functions 
with respect to the corresponding variables and ho is the 
field in absence of a defect, Eq. (11). Evaluation of the 



2D FT of this equation is outlined in Appendix C: 

K p 2 h b (k, z) = -^J ^ Kq - k)r ° A > (42) 

where p 2 = \q 2 + k 2 and the vector A is given by 

A x = Kx - 1 2fl 0x + (k y - q y )(q x h ) z , 
A y = Ky - 1 2h 0y - ( k x - q x )(q x h ) z , 
A z = K z - q 2 h 0z + (q x - k x )(ih' Qx + q x h 0z ) 

+ (Qv ~ k y)( ih 0y + Qyhoz) ■ (43) 



Since ho = He pz with H given in Eq. (11 1 we have: 
-4, 



2,7'.:/' q)q x<y e Pz <p s (q), P 
A z =2(P-q)q-ke Pz <p s (q). 



An 



(44) 



The solution of Eq. ( 42 ) is obtained as described in 
Appendix B: 



h b (k,z) = Ce pz 



rr 



dq e 1 



P 



Ae l{q ' k)ro . (45) 



The conditions of divh-^ = and of the field continuity 
at z — are analogous to Eqs. ( 28) and ( [29]) - ( 31 ): 



ik-C + pC z = 



4tt 2 



ik x ((p s + ip r ) = h 0x + C X - 



ik y (ip s + ip r ) = h 0y + Cy 
-k(tp r - Lp s ) = hoz + C; 



- k 2 


-k)r 


dq A x e^ q - 


-k)r 


47T 2 q 2 - 


k 2 


dq A v e l( - q - 


-k)r 


47r 2 q 2 - 


k 2 


dq A z e^ q 


-k)r„ 


4^2 q 2 _ 


k 2 



(46) 
(47) 
(48) 
(49) 



Using divft-o = and Eq. ( 46 ) we obtain the part of the 



response field due to the defect: 
2r] 4 



Tp{k,r Q ) 



\ 2 



dq {P - q)q-kip s {q)e l (i- k >" 
Att 2 k(p + k)(P + p) 



(50) 



This expression can be considered as the Green's func- 
tion for the general problem of a 2D defect: 

G(r,r )=i>(r,r ), G(k,r ) = ^(fc,r ) . (51) 



The response potential due to a defect distributed with 
the area density N(ro) is 



V(fc) = / dr Q N(r Q )G(k,r Q ) 



(52) 



In particular, for a plane defect situated at Xq = 0, we 
obtain by integrating this over y$: 



5(p r 



KKp + k) 2tt 



dq x (Q - K)k ■ Q 



K 



P 



V s (Q) ■ (53) 



where n is the linear density of the point defects along the 
line x = 0. This coincides with Eq. (32) of the previous 



section and establishes the relation between the constants 
used: /3 3 = nif . 

Another useful example is that of a uniform slab of a 
width w confined between the planes x = ±w/2. The 



outside potential is obtained by integration of Eq. ( 50 ) 
over Xq between ±w/2 and over j/q from — oo to oo: 



r r 2N V 4 

tt\£ 



dg x (Q-K)k-Q sin ^~ k ^ w 

k(p + k)(K+p) q x -k x 

(54) 



It is instructive to use this example to establish the 
relation between the penetration depth of the "defective" 
slab Ad and the constants we are using. To this end we 
write the expression ( 40 ) for a unit area of the slab cross- 



section represented as N point defects: 

N 



(55) 



where r v are the positions of the point defects. Clearly, 
N = S{r — r v ) so that at the slab 



\ 2 = \l-n 4 N. 



(56) 



To relate the factor f3 in Eq. ( 14 ) for the planar defect to 



the characteristics of the slab, one takes the limit w 



in Eq. ( 54 ) and compares the result with Eq. ( 32 ) for 
planar defects to obtain: 

^ = w{\ 2 - A 2 ) . (57) 

Similar to the slab is the case of a cylinder of a radius R 
with penetration depth A^ immersed in a material having 
the penetration depth Aq: 



8tp r = 



rfNR f dq{P - q)q-ktp s {q) ,h{\q-k\R) 



ttA 2 , 



k(p + k)(P + p) 



k\ 



(58) 



where J\ is the Bessel function of the first order. 

Another situation where the Green's function yields a 
straightforward solution for the outside field is a system 
of periodically arranged point defects. For simplicity we 
consider a square lattice of defects with the unit cell size 
d. The integrand in Eq. (|50|) then contains 



E' 



,i(q-fc)r„ _ 



4tt 2 
rJ 2 " 



^2S(q-k-v), 



(59) 



7 



where r n are positions of the defect lattice and u 
are reciprocal lattice vectors having x, y components 
2iri/d, 2itj/d with integers i,j running from — oo to oo 
(see, e.g., Ref.[5J). This gives 



5<p r 



A 0z 



\l<Pk{p - 



k) ^ P 



q—k-\-v 



(60) 



where A 0z is A z of Eq. ( 44 ) at z = 



V. THIN FILMS 

The problem of a linear defect in a thin film with the 
Pearl length Ao = 2A§/d (d is the film thickness, Ao ^> d) 
is formally simpler than for a bulk with a planar defect 
because there is no need to consider the z dependences 
inside the film. 

Let the superfluid density at the y axis of the film at 
z = differ from the rest of the film. The Pearl length 
than can be written as 



A(x) = A - a 2 5(x) 



(61) 



where the constant a, with the dimension of length, can 
be expressed in terms of the superfluid density enhance- 
ment or suppression. 

One can solve the film problem basically along the lines 
described in detail for the bulk case. To avoid repetitions 
we provide here only the result for the magnetic potential 
due to the defect: 



dq x (k ■ Q)(p s (Q,x ) 



k(l + kA Q ) J _ ca 2tt 1+A Q 



(62) 



where Q — {q x ,k y ). (p s (Q,xa) is the 2D Fourier trans- 
form of the magnetic potential of a source at a distance 
xq from the defect line in the absence of a film; for a 
circular current source of the SQUID susceptometer, the 



potential is given in Eq. ( 33 ) 



VI. DISCUSSION 

We have made an effort in this work to develop a for- 
malism to analyze scanning susceptometry data of su- 
perconductors containing planar defects, such as twin or 
grain boundaries, perpendicular to the sample surface. 

Superfluid density on the twin boundaries may, in some 
cases, be enhanced relative to the bulk. 6 Within our 
scheme this corresponds to the parameter ft > 0. In this 
situation, vortices should be repelled by t he b oundary, 
as observed on twinned Ba^ei-zCo^A^P^ 

In most cases, however, the grain boundaries attract 
vortices, in other words, the superconductivity is sup- 
pressed at the boundaries. Within our scheme this is 
described as ft < 0. The suppression of the superfluid 
density on grain boundaries should be observable with 
scanning SQUID susceptometry, but to our knowledge 
this experiment has not been done. 



The Green's function approach developed in Section IV 
may serve as a basis for studying the penetration depth of 
nonuniform materials, one of the outstanding problems 
in applying scanning susceptometry measurements to the 
local determination of A. 

One of the motivations for the current work is that 
although stripes of enhanced susceptibility associated 
with twin boundaries have been observed using SQUID 
microscopy,^ they have not yet been seen in magnetic 
force microscopy. - The failure to observe stripes using 
MFM is puzzling, and it is hoped that the present calcu- 
lations will provide guidance for future investigations. 
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Appendix A: Identities ( 18 | and ( 19 ) 



The first identity is a particular case of the convolution 
theorem for the Fourier transform of a product: 

//•oo 
drer ik - r S(x)f(r) = / dye^fiO, y) 
J — OO 

= |~ dye'^y J ^e^yf(q) = f(Q) , (Al) 

where Q = (q x ,k y ). Similarly, one transforms: 
FT[8'(x)f(r)\ = [ dre- ik - r S'(x)f(r) 



dq x 



'i ~^-( k x- qx)f{q x ,k y ) . (A2) 

Z7T 



To Fourier transform the first term on the RHS of 
Eq. fHi]), we note that x x curlfao = — ycuvl z hQ + 
z curly ho- We then have: 

FT[5'{x)x x curlft-o] = V FT[5' (x)curl z h a ] + 

~^{k x - q x ){iQ x ho) z 
-oo 27r 

-i~ I ^(kx-qx)[ti 0x -iq x ho z ], (A3) 



where h = h {Q,z) and h' 0x = dh (Q,z)/dz. The FT 
of the last term on the RHS of Eq. ( 16 ) is: 



FT[5(x)V 2 h ] 



^[K(Q)~Q 2 ho(Q)]. (A4) 



Appendix B: Solution of the differential Eq. (26) 



The general solution of h"(z) — p 2 h(z) = D(z) reads, 
see, e.g., Ref.[U 

h = de pz + C 2 e- pz 
e pz f z , e~ pz f z 

The lower integration limits here can be chosen arbitrar- 
ily, but this choice affects the constants Ci,2 which are 
eventually fixed by boundary conditions. 

Since we are dealing with a linear differential equation 



(26) with the RHS as an integral (a sum) of the factors 
e"^ , we can take the solution for a particular K and 
then perform the integration (summation). Henc e, we 
set D = e K< > and evaluate the integrals of Eq. (Bl ): 



h = e pz 



C x - 



Mk-p) 



c 2 + 



2p(K + p) 



K 2 -p 2 



(B2) 



Since h should vanish at z — > — oo, C 2 = —l/2p(K + p) . 
The solution becomes: 



^Kz 



Ce pz 



K 2 -p 2 



Ce pz + 



1,2 



(B3) 



where C is a redefined arbitrary constant. 



Appendix C: Integration over the SQUID loop 

Given the FT of the response field h bz (k, zq), one can 
do the integration in real space over the area S of the 
circular SQUID loop: 



$ b = / drh bz (r,z ) 
Is 



dr 



dk 

4^ 



ik-r 



/i 62 (fe,z )e lfc ' 



and do first the integration over the loop of a radius r 
centered at (a;o,0). To this end, one goes from the vari- 
able r = (x,y) to r' = {x — xq, y) centered at (xq, 0): 



dre lkr = e ik * Xo 



is 
Hence, 



dr'e lk -' 



2irr 



J\{kr)e 



2tt 



dk 



Mfc,z )Ji(fcr)e^°. (CI) 



Appendix D: Fourier transform of Eq. (41 1 

The LHS transforms to 

h b {k,z)(l + \lk 2 )-\lh'i. (Dl) 

The 2D FT of the first term at the RHS is readily shown 
to be 

-V 4 J ^[-q 2 h (q,z) + K{q,z)]e^- k >°. (D2) 



The next term transforms to 



dq 



4tt 2 



y{q x hy - q y h x ) 



+ z{ih' x +q x h z )]{q x -k x )e i ^-V r \ 



(D3) 



Here, the arguments (q, z) of all Fourier components hi 
have been omitted for brevity along with the subscript 
denoting unperturbed fields. The FT of the third term 
on the RHS is 



dq 

4^2 



x{q x h y - q y h x ) 



(M y + qy h z )}(q y ~k y )e^-V r °. 



(D4) 



Appendix E: Interaction of a vortex with a parallel 
defect plane in an infinite sample 

An infinite vortex at r v — (x v , 0) perpendicular to the 
sample surface at z — has only the h z (x, y) component, 
with the FT 



h Qz 



f) e 



0o e 



-ik x x v 



1 + X 2 k 2 



X 2 P 2 



(El) 



In the presence of a planar defect at x = 0, A is given 
in Eq. ( 14 ) and the London equation is that of ( 16 ). The 



total field can be written as h z = h$ z + h bz where Hq z is 
the vortex field unperturbed by the twin boundary given 
in Eq. (El ), and h bz is the boundary perturbation. For a 



weak perturbation by the boundary, we obtain: 

h b - \ 2 V 2 h b = -j3 3 6'(x)-^ - p 3 5(x)V 2 h , (E2) 



where ho is given in Eq. (El ) and the subscript z is omit- 
ted. After FT this gives: 



h b (k) 



0o/3 3 
A* 



dq x k ■ Q 

~2n K 2 p 2 



(E3) 



The London energy per unit length is a sum of mag- 
netic and kinetic contributions: 



8tt£ L = J dr [h 2 + A 2 (curlft.) 2 ] 



(E4) 
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Since h = ho + h\, and A 2 = A 2 , — f3 3 5(x), we obtain the 
interaction energy in linear approximation: 



2hohb + 2AgCurlfao • curlhf, 



dr r 

8^ 

f3 3 5(x)(cmlh ) 2 



(E5) 



We now substitute here the Fourier integrals for h and 



>; 20 

D) 

i 

CD 
!= 

® 15 
C 

o 
o 

CO 10 

<5 

H — ' 

c 



0.2 0.4 0.6 0.8 1 1.2 1.4 

X 



hi,, integrate first over r, and take the real part of the 
result: 

c/>q/3 3 f dkdq. 



8tt 4 A^ 



dkdq x ( 

[k x q x sin k x x sm q^a; 



+k y cos /c^a; cos g^x J . (E6) 



Here, the integrand is dimensionless (Ao is used as a unit 
length) and all integrals are from to oo (the integrand 
is even in k x ,k y , and q x ). This integral can be evaluated 
numerically. Fig. [4] shows the resulting repulsive inter- 
action between the vortex and the planar defect with 
enhanced superfluid density, in agreement with observa- 
tions reported in Ref.'sQ]and|9l 

It is worth noting that the calculated energy £i n t{x) 
diverges at x = 0: setting x = and integrating first 
over q x and k x from to oo, one is left with 



7T 

T 



dky ky 
l + kl 



- (to — tan 1 to) 



(E7) 



FIG. 4. (Color online) The interaction energy per unit length 
between a vortex parallel to a planar defect in an infinite 
sample in units of c/>o/3 3 /87r 4 Ao as a function of x = x v /\o, 
evaluated numerically with the help of Eq. ( E6 1 for the super- 



fluid density enhancement at the twin boundary. 



which diverges as to — > oo. The divergence is an artifact 
of our model, which assumes that the effect of the twin 
plane is weak and keeps only linear terms in the correc- 
tion due to the planar defect. At short distances between 
the vortex and the planar defect, the interaction is not 
weak and the model fails. 
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